import numpy as np
import matplotlib.pyplot as plt

def normal(x):
    for idx, item in enumerate(x):
        if item < -100:
            x[idx] = None
    return x

f = np.array([i for i in range(-300, 300, 1)])/100
msk  = (16 / (np.pi**2) )\
                  * (np.cos(6.2832 * f))**2 / (1 - 16 * f**2)**2

BT = [0.3, 0.5]
gmsk = []

for W in BT:
    H_G = np.exp((np.log(2)/(-2.))*(f/W)**2)
    gmsk.append(normal(10*np.log10(msk*H_G)))

msk = normal(10*np.log10(msk))

plt.figure(figsize=(7, 4), dpi=300)
plt.plot(f, msk, label='MSK', color='green')
plt.plot(f, gmsk[1], label='GMSK BT = 0.5', color='red')
plt.plot(f, gmsk[0], label='GMSK BT = 0.3', color='black')
plt.title("Normalized Power Spectral Density")
plt.xlabel('Frequency Offset / Bit Rate (Hz/bit/s)')
plt.ylabel('Spectral Power Level (dB)')
plt.legend()
plt.grid()
plt.savefig('PSD.png')